function [beta2, T, T_J, F, part, OIR] = SURPartiesRobust_IVGMM(Y, X, Z, mean, RAND, n, k, l)

% First estimation
P=Z*inv(Z'*Z)*Z';   %Projection matrix
beta=inv(X'*P*X)*X'*P*Y;   % Get parameter estimates
U=Y-X*beta;  % Get residuals

clear P;
% Second Estimation

P2=Z*inv(Z'*(kron(eye(n),ones(2,2)).*(U*U'))*Z)*Z'; %Second step robust projection matrix
beta2=inv(X'*P2*X)*X'*P2*Y; %Second step parameter estimates
u2=Y-X*beta2;  % Second step residuals
U2=(reshape(u2',2,n))'; % Reshape second residuals
alpha=inv([X(:,2) X(:,k+2)]'*[X(:,2) X(:,k+2)])*[X(:,2) X(:,k+2)]'*u2; % Get alpha parameters for error term
ETA=u2-[X(:,2) X(:,k+2)]*alpha; % Get exogenous residuals
clear U P2;

% Compute COVAR matrix
COVAR=inv((X'*Z)*(inv(Z'*(kron(eye(n),ones(2,2)).*(u2*u2'))*Z))*(Z'*X));

% Conduct Wald test on beta parameters
[F, T, T_J] = Wald_GMM(beta2, COVAR, k); 

% Estimate partial effects
bpartial=[beta2(2:k); beta2(k+2:2*k)]; % Create beta vector excluding constant
mean_l=[mean;mean;mean];    %create matrices for partial effects of binary variables
mean_h=mean_l;
mean_h(1,4)=1;
mean_h(2,5)=1;
mean_h(3,6)=1;

[part] = partial_IVGMM(beta2, bpartial, mean, mean_l, mean_h, RAND, ETA, alpha, n, k);

% Compute overidentification restrictions test statistic

OIR= (Z'*u2)'*(inv(Z'*(kron(eye(n),ones(2,2)).*(u2*u2'))*Z))*(Z'*u2);
